Finite-temperature order-disorder phase transition in a 
cluster model of decagonal tilings 
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Sec. IIIII we extend this model to 3D by stacking the 2D 
layers on top of each other. An additional coupling be- 
tween adjacent layers is introduced. For this 3D model 
we define in Sec. IIVI an order parameter to distinguish 
between the ordered quasicrystalline phase at low tem- 
peratures and the disordered random tiling phase at high 
temperatures. By means of this order parameter, the 3D 
system can then be investigated for a potential order- 
disorder phase transition, which is done in Sec. using 
Monte Carlo (MC) simulations. 
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In a recent paper [Cluster Model of Decagonal Tilings (to be published in Phys. Rev. B)], we have 
introduced a cluster model for decagonal tilings in two dimensions. This model is now extended 
to three dimensions. Two-dimensional tilings are stacked on top of each other, with a suitable 
coupling between adjacent layers. An energy model with interactions leading to a perfect decagonal 
quasicrystal at low temperatures is studied by Monte Carlo simulations. An order parameter is 
defined, and its dependence on temperature and system size is investigated. Evidence for a finite- 
temperature order-disorder phase transition is presented. The critical exponents of this transition 
are determined by finite-size scaling. 
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I. INTRODUCTION 

It is well knowni that, in two dimensions (2D), perfect 
quasicrystalline order cannot be stable at positive tem- 
perature if the interactions have finite range. At positive 
temperature, a 2D quasicrystal is always in the random 
tiling phase without long-range order, the "transition" 
from the ordered to the disordered state being at T = 0. 
For three-dimensional (3D) quasicrystals, on the other 
hand, this order-disorder phase transition is expected to 
occur at positive temperature^ 

The low-temperature state is also called the locked 
phased because its phason degrees of freedom are frozen 
(locked). The high-temperature state is accordingly 
called unlocked phase. Here, the thermal energy is suffi- 
ciently high to excite the phason degrees of freedom. 

3D axial quasicrystals, in particular decagonal ones, 
can usually be regarded as periodic stackings of 2D lay- 
ers, each of which is quasiperiodic. Geometrically, these 
layers can be described as decorations of quasiperiodic 
tilings like the Penrose tilings. Henlej^ has proposed 
to model axial quasicrystals as stackings of 2D tilings 
with a suitable coupling between adjacent layers, in ad- 
dition to the coupling inside the 2D layers. Such layered 
tiling models built up from 2D Penrose tilings have been 
studied by Jeong and Steinhardt^ who indeed found a 
finite-temperature order-disorder phase transition. 

In a recent paper 4 (hereafter referred to as Paper I) , we 
have introduced a cluster model for 2D decagonal tilings. 
In the present article, we will now extend this model to 
3D stackings. Our aim is to investigate an order-disorder 
phase transition, too, but this time in the framework of 
cluster coverings. 

In Sec. m we first give a brief summary of Paper I, 
explaining the important features of our cluster model. 
This includes the principle of cluster density maximiza- 
tion and the ordering of random structures to perfect 
ones by a coupling between overlapping clusters. In 



II. 



TWO-DIMENSIONAL CLUSTER MODEL 
(BRIEF SUMMARY) 



In Paper I, different versions of overlap rules for clus- 
ters in a 2D covering model of the Penrose pentagon tiling 
(PPT) are discussed. The first version, which has a very 
natural realization in terms of a vertex cluster in the PPT 
(Fig. is called the relaxed rule. If we compare this 
cluster with the well known Gummelt decagon^Sii we 
see that it enforces the correct orientation of the large 
B-overlaps, but not the orientation of the smaller A- 
overlaps (Fig.QJ. 

As an alternative to cluster covering as an ordering 
principle for quasicrystals, the principle of cluster density 
maximization is considered as well. A statistical model 
is built where each random tiling is assigned an energy 





FIG. 1: Vertex cluster, superimposed on the Gummelt 
decagon (left), representative A-overlap (middle) and B- 
overlap (right). This cluster enforces the relaxed overlap rule. 
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FIG. 2: Representative example of a flip move in the PPT. 
(The complete set of possible flip configurations is shown in 
Paper I.) 

which is just the negative of the number of vertex clus- 
ters, thus mimicking the cohesion energy of the clusters. 

The density of clusters can then be maximized in a 
MC simulation by simulated annealing, using flip moves 
like the one shown in Fig. [21 The simulation algorithm 
was a usual Metropolis importance sampling scheme: 8,9 
A proposed flip is accepted with probability p = 1 if it 
decreases the energy (AE < 0), but only with probability 
p = c~ l3AE if the energy is increased (AE > 0). 

The states of maximal cluster density are supertile ran- 
dom PPTs (Fig. with an edge length r 2 times that of 
the underlying tiling (where r = (l + yo)/2 is the golden 
number). Naturally, they contain also cluster overlaps 
which do not satisfy the constraints of Gummelt's perfect 
overlap rule, since the relaxed rule does allow for disori- 
ented A-overlaps because of the missing interior orienta- 
tion of the rhombus (Fig. [TJ. 

In a second step, a coupling between overlapping clus- 
ters is added which penalizes such defects. To keep the 
cluster density constant, the simulations are run at the 
level of the supertiling. Each cluster is represented by 
a vertex plus an arrow indicating the orientation of the 
cluster (Fig. 

As MC moves we can still use hexagon flips like the 




FIG. 3: Structure with maximal cluster density. The cluster 
centers form the vertices of a supertile random PPT (gray 
lines). 




FIG. 4: Supertiling with cluster orientations indicated by ar- 
rows. The forbidden A-overlaps, corresponding to tile edges 
with antiparallel arrows, are marked. 




FIG. 5: Creation (top) and shift of defects (bottom) by single 
flips of the two types: hexagon flip (left) and rhombus flip 
(right). 



ones in Fig. [21 provided that the orientations of neigh- 
boring clusters are updated consistently with the under- 
lying tiling (Fig. [3 left). Additionally, on the level of 
the supertiling there is a new type of flip: The rhom- 
bus flip (Fig. right) only changes the orientations of 
the clusters on the obtuse rhombus corners, but keeps 
the tiling itself fixed. In comparison with the basic flip 
moves in the Penrose rhombus tiling 2 ' 9 (which were used 
in the simulations of Jeong and Stcinhardt 3 ), the hexagon 
flip corresponds to D-typc configurations and the rhom- 
bus flip to Q-type configurations of the Penrose rhombus 
tilingia 

As can be seen from Fig. [SJ it is possible to create, 
annihilate, or shift defects by these flips. The defects 
correspond to tile edges with antiparallel arrows at their 
ends (marked in Fig. [SJ). Hence, we assign a positive 
energy to the creation of a new defect and use this energy 
model in a Metropolis MC scheme. 

By simulated annealing it is shown that the cluster 
coupling model is capable of ordering the supertile ran- 
dom PPTs to perfectly quasiperiodic structures at low 
temperatures. Since, in our 3D simulations, we want to 
study transitions between perfect order and disorder, we 
will use this second kind of model for the layers in our 
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FIG. 6: Flip constraint for the interlayer coupling. The flip 
of the interior vertex of the middle-layer hexagon violates the 
congruence of the layers and hence costs energy. 
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FIG. 7: Creation (left) and shift of interlayer defects (right). 
The arrows represent here schematically the orientations of 
hexagons or rhombi, respectively, in the layers. 



stacking of coverings. 



III. EXTENSION TO THREE DIMENSIONS 

We will now consider 3D stackings of our 2D covering 
model. The intralayer interactions are those discussed 
above, which favor Gummelt's overlap rule, and hence 
produce perfectly ordered structures inside a layer at low 
temperatures. In addition, a new coupling between the 
layers has to be introduced. This interlayer coupling is 
chosen such that it prefers adjacent layers to be congru- 
ent. Then the ground state consists of identical, per- 
fectly orderd layers. In the high-temperature state, the 
individual layers are random tilings which need not to be 
congruent, so that there is also disorder in the stacking 
direction. 

Following an idea of Henley^ and analogously to Jeong 
and Stcinhardtj^ we formulate a flip constraint for the 
interlayer coupling. The vertex inside a hexagon in a cer- 
tain layer can be flipped only if there is also a hexagon 
with coinciding boundary (ignoring the interior vertex) 
in the adjacent layers above and below (Fig.[SJ. An anal- 
ogous rule is imposed for rhombus flips: There must be 
a coinciding rhombus above and below in order that the 
middle rhombus can be nipped. These constraints main- 
tain a minimal correlation between adjacent layers. 

As we want to favor congruent layers in the ground 
state, energies are assigned to any mismatch between ad- 
jacent layers. A flip can then either create or remove two 
mismatches, or it can shift a mismatch up- or downwards 
(Fig.0. 



IV. DEFINITION OF AN ORDER PARAMETER 

In order to analyse the transition from the disordered 
phase at high temperatures to the ordered phase at low 
temperatures, a suitable order parameter has to be de- 
fined, which can distinguish between perfect order and 
disorder i&il 

Within one layer, perfect order means the absence of 
disoriented A-overlaps, as discussed in Sec. [HJ In a con- 
tinuous sequence of hexagons and rhombi along a straight 
line (Fig. |SJ, the orientations are then all the same, 
whereas a disoriented A-overlap switches to the oppo- 
site orientation. These lines correspond to the so called 
"worms" or "Ammann lines" in the Penrose rhombus 
tiling, where again the rhombus in the PPT cor- 
responds to Q-type configurations of Penrose rhombi, 
and the hexagon to D-type configurations. Flipping a 
hexagon or a rhombus in a perfect worm creates mis- 
matches and interrupts the sequence of equal orientations 
along the worm line. 

Therefore, if we characterize the orientation of each 
hexagon and rhombus along a straight line by a variable 
Si = ± (Fig. |SJ, we can compare the worm line to a ID 
"spin chain" , assigning "spin up" for one orientation of 
hexagons and rhombi and "spin down" for the other. If 
we sum up all the spin variables Si along a line, they will 
average out to zero in the random tiling phase. In the 
ordered phase, all the spins have the same orientation, so 
that the sum is proportional to the length of the chain. 

We extend this picture to 3D by combining spin chains 
atop each other to "spin sheets" (Fig.^Jl, comparable to a 
2D spin system. In the ordererd phase, all the layers are 
congruent, and thus all the spins within one sheet have 
the same orientation. In the disordered phase, the sum 
over all spins in a sheet will again be zero. However, even 
in the perfectly ordered phase, parallel spin sheets do not 
necessarily have the same spin orientations. Therefore, 
absolute values have to be taken for each sheet separately. 
The order parameter ( "magnetization" ) thus becomes 



{s 3 } ies 3 



(1) 



where n is the number of spins (hexagons and rhombi) 
considered. The inner sum, of which the absolute value 
is taken, runs over all hexagons and rhombi (spins Si) 




FIG. 8: Worm line consisting of hexagons and rhombi. Each 
hexagon and rhombus is assigned a "spin variable" + or — . 
The presence of defects leads to an alternation of the value of 
this variable along the line. 
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21) quasiperiodic plane 

FIG. 9: "Spin chains" one atop the other are combined to a 
"spin sheet". 



lying in the jth spin sheet Sj. Afterwards, the contri- 
butions of all parallel spin sheets Sj are added. Such an 
order parameter is defined separately for each of the five 
directions in the tiling. As we have to use periodic ap- 
proximates with periodic boundary conditions, some spin 
sheets may wrap around the torus several times, which 
has to be taken into account when determining the index 
of the spin sheet on which a given spin is located. 

The value of this order parameter is one in the perfectly 
ordered case, since then all the spins in a sheet have the 
same orientation. In the totally disordered phase, the 
spins average out to zero already along the chains, so 
that the order parameter is zero. The magnetization M 
therefore provides a well defined and suitable order pa- 
rameter for the detection of perfect quasiperiodic order. 

In addition to the magnetization, we also define the 
corresponding "susceptibility" x, which measures the 
fluctuations of the order parameter. In analogy to the 
magnetic susceptibility, we set 



X = n/3 ((M 2 > - (M) 2 ) 
where [3 is the inverse temperature. 



(2) 



V. ORDER-DISORDER PHASE TRANSITION 

The behavior of the order parameter has been investi- 
gated by MC simulations at different temperatures, us- 
ing a series of stacked periodic approximants. In order to 
work out real 3D effects, the number of layers is propor- 
tional to the linear dimension of the approximant. We 
used approximants of order fk+i/fk m the x-direction, 
and of order fk/fk-i in the y-direction. which have the 
least number of defects^ fk is the fcth Fibonacci number 
(fa = 0, fx = 1, fk+i = fk + fk-i)- The linear dimension 
is then of the order of fk-i, so that we took fk-i layers 
in the stacking direction. This resulted in systems with 
total vertex numbers N = 141 (fc = 5), 615 (6), 2576 (7), 
10959 (8), and 46347 (9). The run lengths were of the 
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FIG. 10: System size dependence of the sheet magnetization 
vs. inverse temperature. 
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FIG. 11: System size dependence of the susceptibility vs. in- 
verse temperature. 



order of some 100,000 MC sweeps at each temperature, 
where perfoming one MC sweep means to choose TV-times 
a vertex randomly. The correlation time was found to be 
short compared to the run length, although no precise 
measurement has been made. 

A Metropolis importance sampling schem o 3 i 8 i 9 is used 
for the simulations, similar to the one used in Paper I 
for the 2D simulations. For each proposed MC move, the 
total energy change AE due to intralayer and interlayer 
couplings is computed. The move is then accepted with 
probability p = 1 if the energy is decreased and with 
probability p — e~^ AE if the energy is increased. The 
results for the magnetization and the susceptibility are 
shown in Figs. ITUl and llll respectively, for one represen- 
tative direction in the tiling. 

Fig. illustrates how the magnetization curves con- 
verge to M = 1 at zero temperature (f3 — > oo). Further- 
more, the magnetization at a fixed temperature above 
the critical point ((3 < (3 C , with (3 C « 1.2) decreases with 
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FIG. 12: Finite-size scaling plots L 8/u M vs. L 1/v t for the 
magnetization with the exponents /3 = 0.08 and v = 1.6. 
The critical exponent f3 and hence the scaling function M are 
only defined below the critical temperature, i. e., for the data 
collapse only the range t < has to be considered. 
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FIG. 13: Finite-size scaling plots L" 7/ "x vs. L 1/v t for the 
susceptibility with the exponents 7 = 1.7 and v = 1.6. 



increasing system size. The magnetization curves should 
tend to M — at temperatures above the critical point, 
but they do so only very slowly. This can be understood 
as follows. In the random tiling phase, many spin sheets 
contain only very few spins, so that on those sheets the 
spins do not completely average out to zero. If for the 
magnetization only sheets are taken into account which 
contain a number of spins sufficiently high for a good 
statistics within the sheet, the value of M can be de- 
creased considerably. At (3 = 0, M decreases for the sys- 
tem with 10959 vertices from M w 0.4 down to M w 0.1, 
and even further down to M w 0.08 for a larger system 
with about 75000 vertices. 

The magnitude of the susceptibility, shown in Fig. 1111 
seems to diverge close to the transition temperature with 
increasing system size, which is another evidence for a 



phase transition. The maximum of the susceptibility 
yields a critical temperature of (3 C ~ 1.2 in the limit of 
infinite system size. 

Both the behavior of the magnetization and the sus- 
ceptibiliy are consistent with a temperature-driven tran- 
sition between the ordered (locked) and the disordered 
(unlocked) phase. To obtain the critical exponents of this 
phase transition, the method of finite-size scaling 8 has 
been used to fit the measured curves for magnetization 
M L (t) and susceptibility xl (*) to the scaling functions 

M(LV"t) = L p ' v M L {t) (fart<0), (3) 
X (L^t) = L-^ XL (t), (4) 

where L is a characteristic length of the system (we use 
L = TV 1 / 3 ) and t = T/T c — 1 the reduced temperature 
(with T c = 1//3 C ). The critical exponents j3 and 7 are de- 
fined by the behavior of magnetization and susceptiblity 
close to the transition temperature: 



M oc \T-T c f (for T <T C 
X cx \T-T \-f. 



(5) 
(6) 



The critical exponent v, which determines the behavior 
of the correlation length, £ cx |T — T c | _!/ , has not been 
measured explicitly. By fitting the parameters (3, 7, and 
v in Eqs. © and (@J such that all the curves superimpose 
fFigs. rr21and ll3fl . we obtain the critical exponents 



(3 « 0.08. 
7 w 1-7, 
v « 1.6. 



(7) 

(8) 
(9) 



The uncertainty of these values is about 10% for 7 and v. 
In the case of /3, the spread of values giving a reasonable 
fit is even larger. 

We have also measured the energy E and the specific 
heat C = TV" 1 /? ((E 2 ) - (E) 2 ), where VY is the number 
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of vertices and p the inverse temperature. The result is 
shown in Fig. [21 The maximum of the curve is close to 
the transition temperature determined above. It seems 10 

that the magnitude of the specific heat is independent of A v' 1 ■ !;t ; 46347 

the system size (for large enough systems) , which implies g . i J T \\ 

a critical exponent of a « 0, where a is defined by C cx / -.i'jvj, % , 

\T — T c | _Q close to the critical point. This is in agreement I ,f ' J * ' 

with the scaling relation a + 28 + 7 = 2 (Rushbrooke's h v ; m ^ *f f J ?ji t 

law) A /' /* J l*.t\?^ 

* "*.*-* - 



For comparison, the values obtained by Jeong and 4 - ^ / / 1 j 
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Steinhardti for stackings of Penrose rhombus tilings are 
/3 = 0.2, 7 = 1.6, v = 1.6, and a = 0. These values are in 

good agreement with our result (except for the exponent ,v,' ''* .»•••• ** * ' ' $'"1 

8, whose uncertainty is relatively large in our case). . 5 f ¥ I +--■■+. 
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VI. SOME FURTHER REMARKS 

To test the reliability of our results, we have run sim- 
ulations with different ratios of the intralayer and inter- 
layer coupling strength. It is evident that this influences 
the value of the critical temperature and hence leads to 
different temperature scales. Normalizing the temper- 
ature to the particular critical point, i. e., using B/[3 C 
as temperature scale, we obtain the same behavior of 
the magnetization, independent of the relative coupling 
strengths (Fig. I15|) . Only the curves where one of the 
coupling energies is set to zero show a different behav- 
ior, especially when the energy for the cluster coupling is 
zero, labeled by (0,1). In this case, there is no interac- 
tion inside the layers, so the order parameter is minimal, 
i. e., zero for infinite system size (see discussion of Fig.lTUl 
in Sec. 0). On the other hand, in the case (1,0) where 
the coupling energy in stacking direction is zero, there 
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FIG. 15: Behaviour of the sheet magnetization for different ra- 
tios of intralayer and interlayer coupling strength. All curves 
are for the same approximant. The case of equal coupling 
strengths is labeled by (1,1), the case where the intralayer 
coupling is twice the interlayer coupling is labeled by (2,1), 
and so on. 



FIG. 16: System size dependence of the susceptibility vs. in- 
verse temperature in the case of the relaxed rule. 



is still some kind of purely geometrical coupling between 
the layers due to the flip constraint (Fig. [HJ - This yields 
an order parameter larger than zero. However, it does 
not approach M = 1 at high 3 (low temperatures) since 
this geometrical coupling only restricts the number of 
possible flips but does not enforce congruent layers. 

We have also studied a 3D version of the cluster density 
maximization model obeying the relaxed rule (Sec. ITTll . 
In this case, the intralayer energy is just the negative 
of the number of clusters. (For the interaction in stack- 
ing direction we keep the coupling described in Sec. IHIfl . 
Each ground state now consists of supertile random PPTs 
in each layer, which are all congruent. These supertile 
random PPTs are locally ordered, but show disorder at 
larger scales in terms of disoriented A-overlaps which vi- 
olate the perfect overlap rule. The question is whether 
this partial, local order distinguishes the ground state 
sufficiently from the disordered high-temperature state, 
so that a phase transition between the two can take place. 
We cannot answer this question yet, because the curves 
obtained for the magnetization are too "noisy" (although 
the error bars are small), and the behavior of the suscep- 
tibility (Fig. I16fl does not allow to decide whether or not 
it diverges at the "critical point" , at least not for the ac- 
cessible system sizes. It could well be that the dynamics 
is confined to subsets of the phase space which are sep- 
arated by high energy barriers from other such subsets, 
thus breaking the ergodicity of the simulation. 



VII. SUMMARY 

We have presented a cluster model for 3D decagonal 
quasicrystals which shows a continuous phase transition 
at finite temperature from the ordered low-temperature 
state to the disordered high-temperature state. Using 
the sheet magnetization as order parameter and its asso- 
ciated susceptibility, the critical exponents of this transi- 
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TABLE I: Comparison of the critical exponents 





Q 


p 


7 


V 


Cluster model 


0.0 


0.08 


1.7 


1.6 


Tiling model 3 


0.0 


0.2 


1.6 


1.6 


2D Ising model 15 





0.125 


1.75 


1 


3D Ising model±& 


0.11 


0.33 


1.24 


0.63 



tion have been determined. Within the statistical errors, 
these critical exponents are in good agreement with the 
values obtained by Jeong and Steinhardt 3 for stackings 
of Penrose rhombus tilings (for comparison, see Table QJ. 

The notion of spin sheets suggests a certain resem- 
blance of our model to the 2D Ising model. Indeed, there 
is some similarity also in the critical exponents (TableHJ). 



We have to point out, however, that there is an essential 
difference to the 2D Ising model. The interaction between 
the clusters inside the layers is rather different from the 
interaction between the layers in the stacking direction, 
i. e., the interactions in the system are highly anisotropic. 
Furthermore, in the Ising model only next-neighbor in- 
teractions are considered, whereas in our cluster model 
(and also in the Penrose tiling) the coupling between the 
clusters (or the coupling which arises from the match- 
ing rules in the Penrose tiling, respectively) are of longer 
range. Hence, it is not too surprising that the critical 
exponents of the 2D Ising model deviate from the ones 
obtained for our model (and those obtained by Jeong and 
Steinhardt 3 ). Nevertheless, one can say that the critical 
exponents for the 3D decagonal quasicrystals resemble 
those of the 2D Ising model more than those of the 3D 
Ising model. 



* Present address: Fachbereich Physik, Universitat Kon- 
stanz, D-78457 Konstanz, Germany; Electronic address: 
michael.rcichert@uni-konstanz.de 

1 P. A. Kalugin, JETP Lett. 49, 467 (1989). 

2 C. L. Henley, in Quasicrystals: The State of the Art, edited 
by D. P. DiVincenzo and P. J. Steinhardt (World Scientific, 
Singapore, 1991), p. 429. 

3 H.-C. Jeong and P. J. Steinhardt, Phys. Rev. B 48, 9394 
(1993). 

4 M. Reichert and F. Gahler, Cluster Model of Decagonal 
Tilings (to be published in Phys. Rev. B). 

5 P. Gummelt, Geometriae Dedicata 62, 1 (1996). 

6 P. Gummelt and C. Bandt, Mat. Sci. Eng. A 294-296, 250 
(2000). 

7 P. Gummelt, Combination of Cluster Covering Approach 
and Random Tiling Model (preprint, 2002). 

8 M. E. J. Newman and G. T. Barkema, Monte Carlo Meth- 



ods in Statistical Physics (Oxford University Press, New 
York, 1999). 

L.-H. Tang and M. V. Jaric, Phys. Rev. B 41, 4524 (1990). 
N. G. de Bruijn, Nederl. Akad. Wetensch. Proc. Ser. A 84, 
53 (1981). 

T. Dotera and P. J. Steinhardt, Phys. Rev. Lett. 72, 1670 
(1994). 

J. E. S. Socolar, T. C. Lubensky, and P. J. Steinhardt, 
Phys. Rev. B 34, 3345 (1986). 

A. Pavlovitch, Y. Gefen, and M. Kleman, J. Phys. A: 
Math. Gen. 22, 4347 (1989). 

O. Entin-Wohlman, M. Kleman, and A. Pavlovitch, J. 

Phys. France 49, 587 (1988). 

L. Onsager, Phys. Rev. 65, 117 (1944). 

A. M. Ferrenberg and D. P. Landau, Phys. Rev. B 44, 5081 

(1991). 



